function f = CS(theta2,v,x2,p,delta,ind_market) 

% This function calculates the total consumer surplus following Williams
% (1977), and Small and Rosen (1981)

ns=length(v);   

alpha=theta2(1,1); 
sigmap=theta2(1,end); 
alphai=alpha+sigmap*v(1,:); 

mu1=log(expmu(theta2,v,x2,p));  

expui_i=exp(repmat(delta,1,ns)+mu1);

[xx, yy] = ndgrid(ind_market,1:size(expui_i,2));

expui=accumarray([xx(:) yy(:)],expui_i(:));

%alphai=accumarray([xx(:) yy(:)],alphai_m(:),[],@mean);

f=mean(((0.577+(log(1+expui)))./alphai),2);
end